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ABSTRACT 

We compute the angular power spectrum C, from the BATSE 3B catalog of 1122 gamma-ray bursts 
and find no evidence for clustering on any scale. These constraints bridge the entire range from small 
scales (which probe source clustering and burst repetition) to the largest scales (which constrain possible 
anisotropies from the Galactic halo or from nearby cosmological large-scale structures). We develop an 
analysis technique that takes the angular position errors into account. For specific clustering or repeti- 
tion models, strong upper limits can be obtained down to scales / ~ 30, corresponding to a couple of 
degrees on the sky. 

The minimum-variance burst weighting that we employ is visualized graphically as an all-sky map in 
which each burst is smeared out by an amount corresponding to its position uncertainty. We also 
present separate bandpass-filtered sky maps for the quadrupole term and for the multipole ranges 
/ = 3-10 and / = 11-30, so that the fluctuations on different angular scales can be inspected separately 
for visual features such as localized “hot spots” or structures aligned with the Galactic plane. These 
filtered maps reveal no apparent deviations from isotropy. 

Subject headings : gamma rays: bursts methods: statistical 


1. INTRODUCTION 

The BATSE experiment has now observed more than 
1100 gamma-ray bursts. The observed angular distribution 
is isotropic, while the brightness distribution of bursts 
shows a reduced number of faint events. These observations 
favor a cosmological burst origin. The “great debate” on 
the distance scale of cosmic gamma-ray bursts (GRBs) 
(Fishman 1995; Lamb 1995; Paczynski 1995) considered 
two alternatives; cosmological bursts or events that occur 
in an extended Galactic halo (EGH). The old pardigm of 
nearby Galactic neutron stars with a Population I distribu- 
tion perished due to the combined observations of an iso- 
tropic angular distribution of GRBs along with reduced 
source counts at the faint end of the apparent flux distribu- 
tion (Meegan et al. 1992; Briggs et al. 1996). The absence of 
even a weak “ Milky Way ” band in the GRB distribution 
has indeed made it hard to retain the hypothesis that local 
neutron stars provide the underlying source population. 
Some recent reviews of these and related issues are given by 
Briggs (1995), Fishman & Meegan (1995), and Hartmann 
(1995). 

Although no dominant anisotropies on the sky were 
found in the apparent sky distribution of gamma-ray bursts, 
even small effects might contain valuable information about 
the underlying sources. The detection of a small excess of 
events in special directions, such as nearby stars or the 
Andromeda galaxy, could be a unique signature of stellar or 
Galactic halo models, respectively. For example, a small 
asymmetry with respect to the Galactic plane might suggest 
a local disk origin (Hartmann, Greiner, & Briggs 1995). 


Clustering of bursts beyond that expected from random 
alignments might be evidence of actual clustering of the 
sources or of repeated emission from some sources. Obser- 
vation of repetition would seriously call into question the 
viability of those cosmological burst models that invoke 
unique events, such as mergers of neutron star binaries. On 
the other hand, a detection of the small anisotropy induced 
by the Earth motion relative to the cosmic microwave back- 
ground (CMB) (Maoz 1994; Scharf, Jahoda, & Boldt 1995; 
but see Brainerd 1996) would constitute a convincing proof 
of the cosmological origin hypothesis. These various aniso- 
tropies manifest themselves on different angular scales and 
with different magnitudes. Galactic features would be 
expected to cause large-scale distortions, while burst repeti- 
tion would show its effects on the scale that is typical for 
BATSE source localizations (in excess of 1?6 for the 3B 
catalog). In addition, the instrument does not sample the 
sky uniformly so that we expect some distortions due to the 
nonuniform exposure map of BATSE. 

How should we analyze the angular distribution of 
GRBs? Since the basic null hypothesis of isotropy states 
that burst directions are distributed randomly on the sky 
(which is the impression derived from visual inspection of 
GRB catalogs), we seek tests that can efficiently find 
small deviations. First we search for excess of sources 
toward some direction or a concentration toward some 
plane in the sky, i.e., we seek a dipole of quadrupole 
moment. It is perhaps preferable to search for such large- 
scale anisotropies in an unbiased way by not making refer- 
ence to any particular coordinate frame (Hartmann & 
Epstein 1989; Briggs 1993). On the other hand, such 
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coordinate-free methods are not necessarily the most effi- 
cient ones. If a particular anisotropy is expected, then the 
tests should take this information into account to optimize 
the search efficiency. Paczynski (1990) introduced studies of 
the cos 8 and sin 2 b statistics, where b is the Galactic lati- 
tude of the GRB and 8 is the angle between the GRB direc- 
tion and the vector pointing to the Galactic center. It is now 
common practice to apply both the coordinate-free and the 
Galactic methods to the GRB distribution (Briggs 1993; 
Briggs et al. 1996). These dipole and quadrupole measures 
were sufficient to characterize the large-scale angular 
properties of GRBs when sample sizes were a few hundred 
bursts or less. However, the BATSE experiment has now 
observed so many bursts that an extension of these moment 
methods to higher orders is now useful. In this work, we use 
spherical harmonic analysis (SHA) to represent and inter- 
pret the angular distribution of GRBs 

It can be shown (Horack et al. 1993; Briggs et al. 1996) 
that the statistical estimates of low-order multipoles are not 
very sensitive to the angular smearing induced by statistical 
and systematic localization uncertainties. This is not the 
case for higher order multipoles, which probe the angular 
density field on smaller scales. We shall address this ques- 
tion very carefully in this work. Small angular scales may 
reveal important information about the nature of the GRB 
sources, and location accuracy is crucial. If associated with 
galaxies, we expect clustering on very small scales 
(Hartmann & Blumenthal 1989; Lamb & Quashnock 1993). 
If bursts repeat, we expect clustering at 8 = 0 (Quashnock 
& Lamb 1983a). Both effects are diluted by localization 
uncertainties (the point-spread function), and apparent 
power is transferred from small (or zero) angular scales to a 
scale given by the detector response. A traditional tool for 
the analysis of source clustering is the angular two-point 
correlation function, which was applied first to GRBs by 
Hartman & Blumenthal (1989). The severe reduction in 
correlation strength by positional smearing was demon- 
strated by Hartmann, Linder, & Blumenthal (1991). The 
two-point correlation function is closely related to the 
power spectrum (e.g., Peebles 1980) (in the ideal world with 
no measurement errors or shot noise, one would be found 
to be the spherical Fourier transform of the other). 
However, the correlation function and the power spectrum 
complement each other well, since they are affected by noise 
in quite different ways. This makes it worthwhile to estimate 
both from the data, just as has become the practice with 
galaxy surveys. Another method relevant to the study of 
clustering properties is the nearest neighbor (NN) method 
(e.g., Scott & Tout 1989). This mehod, applied first to GRBs 
by Quashnock & Lamb (1993a), probes only angles near the 
scale defined by the mean angular pair separation, which 
decreases with increasing sample size. We do not consider 
NN methods in this work. 

The remainder of this paper is organized as follows. In 
§ 2, we generalize the standard techniques of power spec- 
trum estimation to properly take into account the location 
errors and the sky exposure of the BATSE catalog. In § 3 we 
apply this to the 3B data set, and in § 4 we discuss the 
results. 

2. METHOD 

In this section, we derive the power spectrum estimation 
technique that is employed in our analysis. The first sub- 


section reviews the statistics of point processes on a sphere. 
This is standard material and has been discussed frequently 
in the literature in connection with the problem of estimat- 
ing the angular power spectrum of point sources such as 
galaxies or quasars (Peebles 1973; Hauser & Peebles 1973; 
Peebles 1980); see Tegmark (1995) for a recent review. The 
extra twist, which makes the analysis of the BATSE data 
more challenging, is the presence of position errors. Since 
some bursts are more accurately localized than others, the 
question of how best to weigh the data is somewhat subtle; 
this is the topic of the second subsection. After that, we 
present the explicit expressions for computing the power 
spectrum estimates from a data set, including a simple beam 
function model. 

2.1. Point Processes on a Sphere 

We model the gamma ray burst distribution as a two- 
dimensional stochastic point process n(r) = JL S(r, r ,), which 
is a Poisson process with intensity (average point density 
per steradian) X(r). Here S denotes the Dirac delta function 
on the surface of the unit sphere, and the unit vectors r r 
correspond to the positions of the various bursts. If we had 
detected a nearly infinite number of bursts, then the func- 
tion 2(f) would be known with great accuracy, and the only 
source of errors when computing its power spectrum would 
be cosmic variance. Since in practice we have only a finite 
number of bursts (in our case 1 122), our estimates of X itself 
will be inexact, leading to the additional complication 
known as shot noise. 

A Poisson process satisfies (see, e.g.. Appendix A of 
Feldman, Kaiser, & Peacock 1994) 

<«(f)> P = m , (i) 

<n(i)n(n> P = mm + m nm . ( 2 ) 

Here X is itself a random field, X(r) = w(r)[l -f A (r)], where 
the underlying density fluctuations A are modeled as a 
Gaussian random field. The function fi, which we will refer 
to as the exposure function , is thus the number of bursts per 
steradian expected a priori, not the number density actually 
observed, In other words, n(f) is proportional to the expo- 
sure time in the sky direction r. As customary, we assume 
that the expectation value <A(r)> fl = 0 and that the sta- 
tistical properties of the field A are isotropic, which means 
that if we expand it in spherical harmonics 1 as 

A (r)= Z Z a lm Y lm (r), (4) 

1 = 0 m - -l 

then 

= & ll Q > (5) 

where the coefficients C { are known as the angular power 
spectrum. There are thus two separate random steps 
involved in generating n: first the generation of the smooth 
field A, then the Poissonian distribution of points. To make 

1 Since all our fields are real valued, we will find it convenient to use the 
real-valued versions of the spherical harmonics throughout. These are iden- 
tical to the conventional spherical harmonics Y lm as defined in, for instance, 
Abramowitz & Stegun (1965), except that the complex exponentials e lnup 
are replaced by 2 1/2 sin (m<p), 1 and 2 l/2 cos (mcp) for m < 0, m = 0 and 
M > 0, respectively. With this convention, the standard identities involv- 
ing spherical harmonics remain unchanged except that no complex conju- 
gation is needed. For instance, the orthogonality relation becomes simply 

J YJj)Y r m .(f)da = 8„ .6^. . (3) 
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this distinction clear, we use the notation <> p for expecta- 
tion values involving only the Poisson step (as in eqs. [1] 
and [2]) and write <>* for expectation values involving only 
the Gaussian field A (as in eq. [5]). When we write simply 
<) without a subscript, we mean the expectation value 
involving both steps. For instance, </i(f)> = <<n(r)) p ) ff = 

< m g = m 

Given the field «(f), we wish to estimate the coefficients 
a im . We denote our estimates a lm , and for reasons that will 
soon become clear, we define them as 

a lm = | Y lm (f) ^ dd - 6 I0 S m0 ^4n . (6) 

We now compute the statistical properties of these esti- 
mates. By substituting equation (1) into equation (6), we 
obtain 

<a lm > = f YJfldR - <5, 0 <5 m0 = 0 , (7) 

i.e., the expectation values vanish. Since the expectation 
values of the true coefficients a lm vanish as well, this means 
that our estimates are unbiased. Notice that we chose to 
include the second term in equation (6) simply to cancel the 
bias arising from the monopole term / = m — 0. Using the 
expressions above, we find that the correlation between two 
multipole estimates is 



X |^<A(r)A(r')> 9 + A. 3(f , t) 2/Q dQ' . (8) 

Substituting equation (4) into this, and using the spherical 
harmonic orthogonality relation (3) and equation (5), this 
reduces to 

\ a im a i'm '/ “ OirVmm’ Q + J wil . (9) 

If h is merely a constant, i.e., if the exposure time is the same 
for all parts of the sky, then the orthogonality relation will 
reduce the second term to simply S ir S mm >/h , and the various 
estimates a lm will all be uncorrelated. Since the true expo- 
sure function h for the BATSE 3B data set varies somewhat 
across the sky, a slight correlation will result. 

Of course, we are also interested in estimating the angular 
power spectrum C z . Defining the quantities 

^lm = &lm ~ him » ( 10 ) 

we find that they are unbiased power estimates (in the sense 
that <£ im > = C z ) if we choose our bias correction to be 


power by averaging the C lm : 




i 

21+1 


l 


Z c, m . 


(12) 


Defining b to be the average of the bias corrections b lm> we 
find that b is in fact independent of /: by substituting the 
spherical harmonic addition theorem (16) into equation (11) 
and using the fact that P x ( 1) = 1, we obtain 


b = 


1 

2/4- 1 


X — 


_t_ 

An 



(13) 


i.e., b is just the spherical average of 1 jh. This means that the 
coefficients b lm , which would be slightly cumbersome to 
compute numerically, need never be computed at all, since 
the power estimate is simply the average of the squared 
a Im -coefficients minus b. 


2.2. The Effect of Position Errors 

The discussion in the previous section applies to any 
population of point sources on the celestial sphere, not 
merely gamma-ray bursts. However, analyzing the BATSE 
catalog involves an extra complication that is absent in, for 
instance, galaxy and quasar catalogs : position errors. 

Let us study first the simple case in which the position 
errors are the same for all bursts in the catalog. If the true 
direction to a burst is f, then we model the apparent direc- 
tion F as a random variable whose probability distribution 
depends only on the angle between r and F. Thus, we can 
write the probability distribution as B(r * F) for some func- 
tion B that we will refer to as the beam function. 

Above, we characterized the distribution of the true burst 
positions as a Poisson process with intensity 2(f), where 2 
was in turn a Gaussian random field. From now on, we will 
let the density n(r) = S (r , f t ) refer not to the true burst 

positions but to the apparent positions. It is easy to show 
that this n will also be a Poisson process, but with a different 
intensity function 2. Specifically, the apparent intensity is 
the true one convolved with the beam function, i.e., 

K PP (?) = (B * x, lu jf) = | Bit • rounder . (14) 

Thus, the effect of the position errors is to smooth out sharp 
features in the expected burst density, which as we will see 
limits our ability to measure fluctuations on scales below 
the beamwidth. Let us expand the beam function in Leg- 
endre polynomials as 

00 (2 1 + l\ 

B(r ■ n = Z o (-^ Jb, P it ■ f) . (15) 


If h is constant, then the bias correction becomes simply 
b lm = l/«, independent of / and m. 

It should now be clear why we divided by h in equation 
(6). If we had not divided by the exposure function, then our 
power estimator C lm would not have measured only what 
we wanted it to, i.e., C v Rather, <£ im > would also have 
received contributions from other multipoles C v , with / ^ 
The quantities C lm are thus good estimates of C x for each 
m-value separately. To reduce error bars, we estimate the 


By using the spherical harmonic addition theorem, 

z rjflYjn = ’ < 16 > 

together with the orthogonality relation (3), we can thus 
write the beam function as 

B{t-t)= £ Z B t YJr)Yjn . (17) 

1 = 0 m= -J 

Applying the beam convolution to equation (4) and using 
the orthogonality relation, we thus obtain the spherical 
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GO / 

(B*AXf)=X I a lm B,Y lm (f). (18) 

1=0 m = -/ 

In other words, convolution with B corresponds simply to 
multiplying the multipole coefficient a lm by B t . 

Repeating the analysis of the previous section including 
position errors (replacing A by B * 2), the case in which h is 
constant 2 thus yields the simple result 

= BfC, . (19) 

In practice, some sources are localized more accurately than 
others, and we clearly want to make use of this fact to make 
the most of the data. Suppose that the total population, 
with number density n, consists of a number of sub- 
populations with number densities n, (so that £ n t = «), and 
that all bursts in the /th subpopulation are equally accu- 
rately localized, as specified by a beam function B { . (Since 
their shape depends only on the sky exposure, all the 
functions h t are identical apart from their normalization, so 
h t x «.) Estimating the power spectrum can now be split 
into two steps: 

1. Estimate a lm separately from each population, as 
above, and call the results a \ m ; 

2 . Combine these estimates into one by some weighted 
averaging, 


< = Z W u aL ■ (20) 

Obviously we want the weights W u to be larger for those 
populations i that are better localized. Let us now deter- 
mine which weighting scheme is optimal. The gener- 
alization of equation (19) to multiple populations is readily 
found to be 


<C lm > = <al>-b l = [ZB li W li ) C,, (21) 

where the bias correction is 


W 2 

b, = Z ^ • (22) 

i n, 

How should we choose our weights W ti l First of all, to 
make C lm an unbiased estimate of C h clearly we want to 
normalize the weights so that the expression in parentheses 
in equation (21) equals unity, i.e., so that £,• B u W Vl = 1. 
Second, we want the error bars on our estimate to be as 
small as possible, i.e., we want to minimize the variance of 


2 When n is not constant, the B * (nA) term in addition gives rise to a 
weak mode coupling between the different multipoles. As discussed below, 
the n of the BATSE 3B data set is basically constant except for small dipole 
and quadrupole corrections. This means that a rm will pick up small contri- 
butions from a rm , where \ l' — l \ <2, which is completely irrelevant for this 
analysis. The reason is that it is merely a second-order effect: we are 
investigating whether there is any signal at all apart from the shot noise, 
and this coupling effect would only alter the relative level of the signal by a 
few percent. Thus, the only instance in which the anisotropy of h must be 
taken into account is when computing the noise bias with eq. (11), since an 
error of a few percent in the (much larger) shot noise contribution could be 
of the same order as the weak signal we are trying to detect. 


C { . In the approximation that a lm is Gaussian, 3 we have 
simply V(€f) = 2(af m } 2 , so that we minimize the variance by 
minimizing the expectation value <a 2 m > = C, + b h Since C, 
is independent of our weights, we thus wish to choose W H so 
as to minimize the bias correction b h given the above- 
mentioned normalization constraint ^ B fi W n =\. This 
constrained optimization problem is solved readily by the 
method of Lagrange multipliers, and the solution is 

W li = b l h i B li , (23) 

where the minimal bias correction is 

b, = (Z B fij < 24 ) 

In summary, we have found our best multipole estimate to 
be 



(For brevity, we omit the trivial monopole correction eq. 
[7] here and below, since we are never interested in comput- 
ing the monopole anyway.) 


2.3. Power Spectrum Estimation in Practice 

For any given data set, the density field n { is just a sum of 
delta functions, one for each burst, so equation (25) reduces 
to 


a 


lm 


Nbj 

4tt 


,Vj 


Z B U z 


w 

n(fj) ’ 


(26) 


where N- t denotes the number of bursts in the /th sub- 
population. We can simplify this expression further by a 
mere change of notation. We let the index k refer to sums 
over the entire burst sample (k = 1 , . . N), and from here 
on, we simply let B lk denote the beam factor corresponding 
to the subpopulation that the kth burst belongs to. Then 




Nf 

4n 


(27) 

where we have defined the effective number of bursts at a 
given multipole as Nf f - { B? k . With this same conven- 

tion, replacing the double sum over subpopulations and 
their members by a single sum over all bursts, equation (25) 
simplifies to 


— 


N 

Nf 1 


Z B*~, 


r, m {h) 


n(?k) 


(28) 


Thus, we have eliminated the need to keep track of 


3 The Gaussian approximation is good when the number of bursts is 
large, by the central limit theorem. It should be emphasized that even 
under circumstances in which this approximation is poor, our C, will be a 
good estimate of the power spectrum: it will simply have slightly larger 
error bars than it would with optimal weighting. 
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subpopulations altogether, 4 and expressed our multipole 
estimates directly in terms of the observed quantities. 

Repeating the analysis for an arbitrary exposure function 
fi, equation (27) becomes generalized to 


In the limit o k 1 (valid for all bursts in the sample as 
shown in Fig. 2), we have to a good approximation that 

(34) 


N 1 CdQ 
' l “ Nf An J n(f) ‘ 


(29) 


We estimate the C x by averaging over m-values as before, 
i.e., 


'' \2l + 1 


~b, 


(30) 


The position uncertainties A 6 quoted in the BATSE 3B 
catalog are defined as the radius of the 1 a circle, i.e., of the 
circle that contains erf [2 -1/2 ] % 68% of the probability. 
Thus, in the limit o k 1, the conversion between A 0 and o 
is 



In the above-mentioned Gaussian approximation, the C lm 
of equation (10) are almost independent with variance 
V(C lm ) = K(af m ) = 2 (a/ m ) 2 , since the b lm are mere constants. 
Hence, the 1 a error bar is 


Note that the values of A 8 quoted in the BATSE 3B catalog 
do not include the systematic error contribution of 1?6, 
which is to be added to the quoted values in quadrature. 
This yields the distribution shown in Figure 2. 


AC, 


1 

2 / + 1 



21 + 1 


1/2 

(C, + b t ) . 


(31) 

Thus, as l increases, the error bars will typically decrease 
first due to the growing number of independent m - modes 
and then gradually start increasing again around the scale 
corresponding to the position errors as Nf { eventually 
approaches zero, making b t explode. 


2.4. The Beam Function 

We model the BATSE beam function as a Fisher function 
(Fisher, Lewis, & Embleton 1987): 




exp (q ft 2 r • r) 
Anal sinh (a k 2 ) 


(32) 


characterized by a location error o k . This is often con- 
sidered by mathematicians to be the spherical version of the 
Gaussian distribution, and it reduces to 


B k (cos 0) « e x p 1 E .V; 2na 2 k , (33) 

when c k 1 radian ss 60°. The Fisher function has the 
advantage that it is correctly normalized (its integral over 
the sphere is unity) for arbitrarily large angles o k , which is 
not the case for the plane Gaussian of equation (33). It 
should be emphasized that although the BATSE location 
error distribution has usually been modeled as a Gaussian 
distribution, it is currently not well enough known that one 
particular distribution is preferred over another, so the 
choice is merely one of convenience. 


4 The Gaussian assumption that we used for computing error bars was 
strictly valid only when 1 for each subpopulation. However, since the 
BATSE 3B distribution of position errors forms a smooth continuum, we 
expect the error bars derived from the Gaussian approximation to remain 
accurate anyway, as long as Nf f $> 1, and this is indeed confirmed numeri- 
cally by Monte Carlo simulations. We generated 1000 mock BATSE 3B 
catalogs with no clustering and analyzed them with the same software as 
the real data, extracting the multipoles / ^ 40. To within the Monte Carlo 
errors (a relative error of order 1000“ 1/2 % 3%), the actual error bars were 
identical to those expected analytically when making the Gaussian approx- 
imation. 


3. RESULTS 

We have used the improved BATSE positions of the 3B 
catalog (Meegan et al. 1995b, 1995c) to expand the angular 
distribution of GRBs in terms of spherical harmonics. The 
3B catalog contains 1122 bursts with known best-fit posi- 
tions (shown in Fig. la) and their statistical uncertainties. In 
addition to statistical shifts, we must also include (in 
quadrature) a 1?6 systematic uncertainty. This value is sig- 
nificantly lower than the 4° of earlier catalogs, and it allows 
us to extend spherical harmonic analysis to / — 50-60 
before localization uncertainties completely wash out any 
possible intrinsic angular power in the GRB sky map. The 
distribution of the actual statistical errors is shown in 
Figure 2. 

Because the sky exposure of BATSE is not uniform 
(Fishman et al. 1994; Meegan et al. 1995b), artificial 
moments are induced (e.g., Briggs et al. 1996). The BATSE 
experiment does not exclude any area of the sky, but due to 
blocking by the Earth and detector gaps during passages of 
the South Atlantic Anomaly (SAA), some positions on the 
sky have a reduced probability for burst detection. The 
associated exposure map is thus best described as a semi- 
transparent mask. While the exposure corrections are not 
as severe as those encountered in galaxy surveys, it should 
and can be included in the analysis. We shall discuss the 
effect of uneven sampling in the next section. 

3.1. The Exposure Function 

Because of problems due to the loss of the spacecraft tape 
recorders, the absolute efficiency has not been determined 
since the release of the IB data set. However, the shape of 
the exposure function h is essentially independent of time, 
and since the shape is all that matters for the present 
analysis, we employ the IB estimate (Fishman et al. 1994). 
This function h depends on declination only and is indepen- 
dent of right ascension. This means that in equatorial coor- 
dinates, the multipole coefficients h lm vanish except when 
m = 0. The dominant deviation from uniformity is a quad- 
ruple (h 20 /h 0 o « 8.8%) depletion of bursts near the 
equator due to the shadowing of the sky by the Earth. The 
second largest anisotropy is a dipole moment (h 10 /h 00 « 
4.5%) toward the Earth’s north pole, due to the South 
Atlantic Anomaly, which requires disabling triggers. Com- 
pared to the shot noise, the higher multipoles (/ > 3) are 
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Fig. 1. — The BATSE 3B data set and the smoothed burst map. The measured locations of the BATSE 3B sample of 1 122 gamma-ray bursts are shown in 
Hammer-Aitoff projection in galactic coordinates {top) and with each burst smeared out by an amount corresponding to the uncertainty in its position 
{bottom). 


negligible (a l0 /a 00 < 1%), but for completeness, they have 
nonetheless been included in our analysis. 

3.2. The Power Spectrum 

The power spectrum C { extracted from the BATSE 3B 
data set is shown in Figure 3, and as can be seen, there is no 
evidence of deviations from isotropy on any angular scale. 
What is plotted is, of course, the difference between two 
positive quantities, the power in the data minus the bias 
correction, according to equation (30), which is why some 
(unphysical) negative estimates occur. Thus, if the gamma- 
ray bursts are in fact completely uncorrelated, we would 
expect the points in Figure 3 to be scattered symmetrically 
around zero, with roughly equal numbers above and below 
the horizontal axis, and about 68% within the shaded 
region. Since all power is by definition positive, the presence 
of any type of correlation would shift the distribution 
upward, leading to a positive excess. 

In Figure 3, we have divided the power spectrum by An to 
make the interpretation of the numbers simpler. A mono- 


pole CJAn — 0.0001 would simply correspond to a fluctua- 
tion of 0.0001 1/2 = 1% in the average burst density. 
Likewise, (C,/4 tt) 1/2 can be interpreted as the density fluc- 
tuation on the angular scale 6 « 60°//. 

Let us comment briefly on this factor of 60 and the 
correspondence between / and 6. From equation (34), we see 
that roughly speaking, a burst probes the multipole / only if 
the factor B lk is of order unity, i.e., if a k l < 1. Here o k is 
measured in radians, so since 1 radian is 1 80 °/tc « 57°, this 
means that only bursters with a location error a < 60 °/l are 
sensitive to the multipole /. 

3.3. The Error Bars 

The size of the error bars (the height of the shaded region) 
in Figure 3 is readily understood from equation (31). For 
/ — 0, we have = N = 1 122, so apart from the factor of 
2 1/2 , the shot noise gives just the familiar Poisson variance 
1 /N. As l increases, the (21 + 1) denominator reduces the 
error bars, since many independent modes are being aver- 
aged. However, as / increases beyond the scale correspond- 
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Fig. 2.—A histogram of the position errors A 6 is shown for the BATSE 
3B sample of 1122 gamma-ray bursts. The 1?6 systematic errors are 
included here, added to the statistical errors in quadrature. 


ing to the typical location errors, the sharp drop in Nf ( 
causes the error bars to increase dramatically. Thus, we 
cannot place strong constraints for / 60 simply because 

there are no bursts that are localized to better than 1?6 



0 10 20 30 40 50 60 

1 

Fig. 3. — The shot-noise-corrected angular power spectrum. Filled 
squares show the multipoles estimated from the BATSE 3B data set with 
minimum-variance burst weighting and shot noise removed (this is why 
unphysical negative values occur). The shaded region shows the 1 a shot 
noise error bars, so if there is no clustering whatsoever, about 68% of the 
squares would be expected to fall within this region, distributed symmetri- 
cally abound zero. Any type of clustering would drive the points upward, 
leading to more points above zero than below. The double-shaded region 
shows that the error bars would be if there were no position errors. 


° 0 20 40 60 

1 

Fig. 4. — The factor by which fluctuations are suppressed by the effect of 
position errors, 2Vf ff /JV, is plotted as a function of multipole /. Our method 
corrects for the smearing by dividing by this suppression factor, which is 
the reason that the error bars in Fig. 3 explode for large l. The suppression 
factor for the real data ( shaded} is compared with the hypothetical situation 
in which all bursts have the same position errors AO, taken to be 1?6 ( solid 
line), 4° (long -dashed line), and 10° (short-dashed line). 

(/ ~ 35), so / > 60 would be more than “2 a” out in the 
Gaussian tail of B h causing the shot noise to explode. 5 

This effect is the reason that the actual error bars become 
so much larger than the “ideal world” error bars (double- 
hatched) that would result if there were no position errors. 
This is also illustrated in Figure 4, where Nf { is plotted as a 
function of L For l = 30, for instance, we are effectively only 
making use of about 25% of all bursts, the remainder being 
too poorly localized to contribute much information about 
the power on this small a scale. Conversly, Figure 3 shows 
also that for / < 5, the location errors have little impact on 
the error bars, confirming the results of Horack et al. (1993) 
and Briggs et al. (1996) for dipole and quadrupole moments. 

Note that Nf* in Figure 4 is far from being Gaussian : for 
small /, it falls off roughly as a Gaussian with A 6 = 10°, but 
for larger /, the tail falls off much more slowly, since most of 
the contribution is coming from the best localized bursts. It 
should also be noted that since the C r coefficients are rota- 
tionally invariant quantities, Figure 3 would look identical 

5 If the true location errors should turn out to be larger than those we 
have assumed here, then the error bars would thus explode at lower /- 
values. For instance, Graziani & Lamb (1995) analyzed the distribution of 
3B locations in comparison to the positions derived from the IPN 3 
network and conclude that the systematic error of the 3B data should be 
about 4° instead of the advertised 1?6. In addition, there may be corre- 
lations in the data that suggest a brightness dependence for the systematic 
error, instead of the constant value suggested by Meegan et al. (1995c). The 
studies by Graziani & Lamb (1995) did not take systematic effects in the 
IPN localization method into account and also do not incorporate the fact 
that some IPN locations are based upon the earlier BATSE 2B locations 
and thus may be biased against the 3B locations. Although we tend to 
agree more with the error budget prescirbed by the BATSE Team, this is 
an unsettled question, and it should be borne in mind that whereas our 
multipole estimates for / 15 (corresponding to 0 > 4°) are virtually inde- 

pendent on this controversy, the estimates for l > 30 should be taken with 
a grain of salt if the reader favors a 4° systematic error. 
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if Galactic rather than equatorial coordinates had been 
used when generating it. 

3.4. Is the Exposure Function Correct ? 

If the estimate of h (Fishman et al. 1994) were incorrect, 
this could introduce artificial signals into our power spec- 
trum. Because of the azimuthal symmetry, this would affect 
only those coefficients a lm that have m — 0. These are 
plotted in Figure 5. Thus, if the bursts are uncorrelated and 
the h estimate is correct, the points should scatter symmetri- 
cally around zero with about 68% of them in the shaded 
region, which appears to be the case. Figure 5 thus provides 
reassuring evidence that h has been correctly modeled. To 
indicate the sensitivity of this analysis, the figure shows also 
the dipole and quadrupole that would be expected if we had 
failed to correct for the above-mentioned Earth-shadow 
quadrupole and the South Atlantic Anomaly. Since the 
quadrupole correction was about 9%, this shows that 
uncertainties in the modeling of the higher multipoles of 
which are typically at least an order of magnitude smaller, 
will not be important compared to the errors 

caused by shot noise. 

3.5. The Minimum-Variance- Weight Burst Map 
Using equation (17) and the orthogonality relation (3), 
whe can rewrite equation (28) as 

= Y lm (f)x(f)dn , (36) 

where we have defined x, the smoothed burst map , as 


N 


I 

k - i 


B\r • f k ) 
n(r) 


(37) 


Thus, we see that the minimum-variance method we derived 
above has a very simple interpretation: apart from the 



Fig. 5. — The multopole coefficients with m = 0. The multipole coeffi- 
cients a l0 in equatorial coordinates, corresponding to fluctuations indepen- 
dent of right ascension, are shown (filled squares ) together with the 1 a 
region expected from shot noise alone. For any isotropic fluctuations, the 
distribution should be symmetric around zero. The triangle and the star 
show the effect that the South Atlantic Anomaly and Earth shadowing 
would have if they were not taken into account in h. 


overall weighting factor N/Nf f , our optimal estimates of 
the multipoles a lm were just the spherical harmonic coeffi- 
cients of a map in which each burst is smeared out by its own 
beam function , and corrected for the uneven sky exposure. 
This map is shown in Figure 1 b and Figure 6 (Plate 1) 
(upper left). A comparison of this map with that using 
earlier BATSE data (Hartmann et al. 1994) shows the tre- 
mendous improvements due to the reduction of systematic 
position uncertainties from 4° to 1:6 and the increase in 
sample size. 

It is quite useful for inspecting the data set visually, since 
in a sense it displays only the information that is really 
present in the data and not more. It does not mislead the 
eye by exaggerating the accuracy to which the burst loca- 
tions are known, enabling those bursts that are well local- 
ized to visually stand out against the background. 

3.6. Bandpass-Filtered Maps 

Although the angular power spectrum C, provides a 
useful measure of the amount of clustering on different 
angular scales, it should be borne in mind that it does not 
contain any information about the relative phases of the 
different multipoles a lm . The same can be said about the 
correlation function, a useful statistical quantity that has 
been estimated elsewhere (Meegan et al. 1995b, 1995c; 
Blumenthal 1995). The loss of phase information means 
that although the power spectrum may tell us that there is 
extra power on some scale, it does not tell us anything 
about where in the sky this power is coming from; we may, 
for instance, be interested in knowing if there are any signals 
localized near the Large Magellanic Cloud or aligned with 
the Galactic plane. Fortunately, this type of information 
(which can be seen as complementary to that provided by 
the power spectrum) is easy to extract with the formalism 
developed above. We define x,(r), the multipole map corre- 
sponding to multipole /, as the sky map 

Xtf) = X a, m Y lm (?) , (38) 

m — —l 

where the estimated spherical harmonic coefficients a lm are 
those defined by equation (28). Similarly, we define the 
bandpass-filtered map corresponding to a multipole range 
/i < / < / 2 a s the sum of the multipole maps for the different 
/ values in the range. Figure 6 shows the filtered maps 
corresponding to / = 2 (the quadrupole), / = 3-10, and 
1= 11-30, respectively, and the reader is encouraged to 
scrutinize these images in search of any features that are 
spatially localized or aligned with the Galactic plane, both 
of which would provide evidence against isotropy. The 
quadrupole, for instance, is neither aligned with the Galac- 
tic plane nor with the equator of Earth, and as is seen in 
Figure 3, its amplitude is of the order that is expected from 
mere shot noise fluctuations. 

Using the orthogonality relation, we see that apart from 
the shot noise correction and a proportionality constant, 
our multipole estimate C { is just the integral of the square of 
the corresponding multipole map, { xf dQ. It is in this sense 
that the filtered maps allow us to see where the power (the 
fluctuations) is coming from. Also, apart from normal- 
ization issues (for instance, the density modulation in h is 
eliminated in the filtered maps), the smoothed burst map in 
Figure 1 h is just an average of all the multipole maps, 
weighted by inverse noise level. Thus, we can think of the 
filtered maps roughly as a decomposition of the smoothed 
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burst map into its different frequency components, into its 
contributions from different angular scales. 

4. DISCUSSION 

Much of the current debate on the origin of GRBs rests 
upon a careful analysis of their angular and brightness dis- 
tribution. Without established counterparts or other burst 
properties that could be used to estimate distances, we do 
not even know their distance scale, which in turn leaves 
burst energetics undetermined. Building models is a chal- 
lenge under such conditions. One of the most important 
pieces of information that we can obtain is the angular 
distribution of GRBs. Deviations from isotropy on some 
angular scale for some or all bursts could provide crucial 
hints to the distance scale. The lack of large anisotropies 
makes it very hard to retain traditional models of neutron 
stars in the Galactic disk. But even models that invoke a 
very extended halo do predict small anisotropies that 
should emerge eventually from the data. And while cosmo- 
logical models generically result in isotropic distributions, 
they too may have tell-tale deviations. We may consider the 
small deviations due to the Earth’s motion with respect to 
the CMB, or the granularity due to local superstructures in 
the cosmic mass distribution. In addition, the well-known 
angular correlations of many cosmological objects or clus- 
tering that would result from burst recurrences would lead 
to some deviations from isotropy. The distribution of burst 
positions on the sky could be the primary source of infor- 
mation leading to an understanding of the burster distance 
scale, and perhaps their nature as well. 

The crucial objective of our study is thus an advanced 
analysis of GRB positions. There are two significant steps in 
this field : 

1. Providing accurate locations for all bursts; 

2. Analyzing this position information with appropriate 
statistical tools. 

The BATSE Team has made great progress in the first area, 
now providing location accuracies down to about 2° for 
many bursts and about 5 C for the average burst (Meegan et 
al. 1995c). The reduction of systematic uncertainties is 
essential for studies of small-scale anisotropies, but it also 
contributes to better estimates of more global patterns that 
may be present in the data. The smearing of burst positions, 
unavoidable from the instrumental point of view, must be 
included in the data analysis. Additional features that must 
be accounted for are temporal and angular gaps in the 
observations. Here we do not consider possible structure in 
burst arrival times, but we study exclusively their arrival 
directions. The exposure function of BATSE must be and 
has been included in this work. 

The remaining question is about selecting appropriate 
tools. This depends somewhat on the question we wish to 
address. Global anisotropies present in many Galactic burst 
models can be studied through low-order multipole expan- 
sion, e.g., dipole-quadrupole statistics, while clustering is 
generally approached with angular correlation functions or 
nearest neighbor distributions. Because of the larger data- 
base and the superior position accuracy of the 3B data 
studied here, we are actually able to bridge these two dis- 
tinct approaches by extending dipole and quadrupole 
analysis of the angular distribution of GRBs to higher order 
multipoles. The technique is the well known spherical har- 
monic analysis (SHA), i.e., expansion of the burst distribu- 


tion in terms of spherical harmonics, y,n,- As discussed 
above, the angular scale probed by a given harmonic is 
approximately 60°//, so that the expansion must be carried 
out to /-values in excess of 30 if we wish to reach the smear- 
ing scale of the current BATSE configuration. 

4.1. Limits on Repetition 

If some fraction of the observed GRBs are repeat events, 
the sky distribution should show angular concentrations on 
small scales (roughly given by the beam smearing of the 
instrument). Evidence for burst recurrence was found in the 
IB data (Quashnock & Lamb 1993a), but subsequent 2B 
data did not confirm this result (e.g., Meegan et al. 1995a). 

The 3B data set is greatly improved over the 2B data in 
its ability to test the repeater hypothesis for the following 
reasons: 

1. The systematic position uncertainty has been reduced 
from 4° to 1:6, and 

2. In addition to the overall exposure time being 
increased by about a year, the post-2B portion of the 3B 
catalog has a greater fractional exposure (live time), which is 
important for repeater models in which the bursting phase 
of sources is less than the BATSE lifetime. 

Burst recurrence is expected to generate excess correlations 
at 6 = 0, which corresponds to excess power at all multi- 
poles. 6 Our study does show some modes with deviations 
around the 2 <r level, but this is by no means a significant 
excess of power because only about the expected number of 
points deviate by about 2 a and the points are generally 
scattered with 1 a of no power. The data are consistent with 
the hypothesis of no recurrences. The althernative hypothe- 
sis tested most frequently in the literature are repeater 
models in which a fraction / of all observed bursts can be 
labeled as repeaters that are observed to burst v times each. 
Tegmark et al. (1996, hereafter T96) employed an SHA- 
based technique to test this two-parameter family of models 
against the BATSE 3B data, and they find that all models 
with (v — 1 )/> 0.05 are ruled out at 99% confidence, as 
compared to the best previous 99% limit (v — 1 )/> 0.27 
(Meegan et al. 1995c). Thus, even a cluster of six events from 
a single source would have caused excess power above that 
present in the 3B catalog. In other words, the multipole 
information that our SHA extracts from the data, as plotted 
in Figure 3, translates into sharp quantitative limits on 
repetition. 

It is conceivable that bursts repeat once or more often on 
a timescale of months and become dormant afterward for a 
much longer period. In that situation, accumulation of 
bursts into a growing sample would dilute the repeater 


6 From the additional theorem (16), one obtains the well-known result 
that 

<A(r)A(r')> = X ‘ P '> C ‘ ■ (39) 

i.e., the C r values are basically the spherical Fourier coefficients of the 
angular correlation function. Correlations only at zero anglular separation 
(before position errors are added in) correspond to the correlation function 
being a Dirac delta function. Just as the regular Fourier transform of S(x) is 
a constant, the power spectrum corresponding to repeaters is C, constant, 
independent of l 



PLATE 1 



Pi G 6. All-sky maps of the 1122 gamma-ray bursts in the BATSE 3B data set, in Hammer-Aitoff projection. In the smoothed burst map ( top left), each burst is smeared out by an amount 

corresponding to the uncertainty in its position. The other maps are bandpass filtered, showing the quadrupole (/ = 2, top right), intermediate scales (3 < / < 10, bottom left), and small scales (11 < / < 30, 
bottom right). 

Tegmark et al. (see 468, 221) 
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signal. When the 3B set is divided into four sets of roughly 
equal number of bursts (not equal time), the correlation 
function shows some small-angle excess at the approx- 
imately 1-2 o level in all subsets (Meegan et al. 1995b; 
Blumenthal 1995). Adding these correlation functions 
together generates a noticeable, but still not highly signifi- 
cant, excess of burst pairs near about 5\ Our corresponding 
SHA analysis for the four subsets (Fig. 7) reveals this effect 
also, but it is evident from this figure that the significance of 
this increase is marginal at best. In other words, SHA yields 
results that are consistent with those obtained by corre- 
lation function analysis. This emphasizes the fact that the 
SHA method now bridges the range of power estimators 
previously employed in GRB studies. 

4.2. Limits on Large-Scale Clustering 

Angular power spectra also constrain burst models that 
trace the large-scale structure of the universe. If GRBs trace 
the galaxy distribution (as neutron star binary mergers 
would), we expect to find angular correlations similar to 
those observed for galaxies or clusters of galaxies 
(Hartmann & Blumenthal 1989; Lamb & Quashnock 1993). 
However, if BATSE samples to a redshift of order unity 
(assuming standard cosmology and standard candles for 
bursts), the sparse sampling of the galaxy density (specific 
rate of bursts inside or near galaxies is ~~ 10“ 6 yr _ l ) and the 
imperfect angular resolution reduces the expected strength 
of the clustering signal. With increasing sample size, it will 
be possible to apply brightness selections, while retaining 
good angular resolution. So far, only the dipole and quad- 
ruple term have been investigated as a function of appar- 
ent source brightness and interpreted in the context of 



Fig. 7. — Angular power spectrum when grouping into quarters. Filled 
squares in the bottom plot show the multipoles estimated as in Fig. 3, 
except that the data have been split into four sequential quarters according 
to when the burst occurred. Each square shows the average of the four 
estimates that the multipole, and the corresponding error bars are seen to 
be twice as large as in Fig. 3, shown above for comparison. The slight 
apparent excess of power in the bottom figure is consistent with the corre- 
lation function analysis of Meegan et al. (1995b„ 1995c), which finds weak 
clustering when the data is time binned into four quarters. 


Galactic GRB models (Quashnock & Lamb 1993b; Gure- 
vich et al. 1993, 1994) and cosmological GRB models 
(Hartmann, Briggs, & Mannheim 1995). 

Figure 3 shows that the BATSE 3B data are consistent 
with the null hypothesis of no clustering at all. Translating 
this multipole information into quantitative limits of course 
requires a detailed specification of a clustering model. For 
the above-mentioned repeater case, this was done by defin- 
ing a two-parameter family of repeater models and placing 
lower limits on the parameter values (T96). Although any 
type of clustering would produce excess power on some 
scale, addressing the issue of how high above the Poisson 
noise this excess would be, and thus the issue of whether 
such clustering is ruled out by BATSE or not, would require 
detailed modeling of the clustering pattern in question. 
Since SHA tightened the previous repeater limits by a factor 
of 5, it appears quite promising to carry out such modeling 
of other clustering scenarios as well. 

If a model predicts clustering that is anisotropic (for 
instance, an excess in the supergalactic plane), then it is 
likely that even stronger constraints can be obtained by 
making use of this spatial template information, since our 
averaging over m -\ alues can “dilute” an anisotropic power 
contribution by up to a factor of (21 4- 1). In other words, if 
one knows exactly what kind of clustering one is looking 
for, one should adopt a data analysis technique that incorp- 
orates this information. In constrast, since the angular 
power spectrum involves no preferred directions in the sky, 
SHA is a useful general-purpose diagnostic tool which will 
detect clustering on any scale (if it is strong enough to stand 
out over the Poisson noise), without requiring any prior 
knowledge as to what kind of clustering one is looking for. 

4.3. Outlook 

We conclude that multipole expansion of the projected 
distribution of GRBs does not show evidence for clustering 
on any angular scale. This argues against the recurrence of a 
substantial fraction of burst sources (at 99% confidence, not 
more than 5% of the BATSE 3B bursts can be due to 
repeating sources [T96]) and against any source population 
with intrinsically strong anisotropies resulting from an 
intrinsically special position of the observer. The remark- 
able degree of isotropy of GRBs constrains severely any 
burst model that invokes traditional geometric features of 
the Milky Way (disk, bulge, or halo). If one wishes to retain 
the Galactic origin hypothesis by introducing very extended 
halo distributions, it seems that these populations cannot 
contribute significantly to the dynamics of the Galaxy 
(those that do are all known to be highly anisotropic) but 
must constitute a trace population. Whether high-velocity 
neutron stars injected into the Galactic halo can indeed 
provide the necessary isotropy remains to be determined, 
and model builders should verify that the proposed spatial 
distributions indeed generate essentially zero angular power 
on all scales, as our analysis suggests. The currently fashion- 
able paradigm of cosmological bursts now passes this test, 
but eventually some deviations from isotropy are expected, 
and spherical harmonic analysis is a tool well suited to 
detect such deviations. 
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